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Dynamical quantum-cluster approaches, such as different cluster extensions of the dynamical 
mean-field theory (cluster DMFT) or the variational cluster approximation (VCA), combined with 
efficient cluster solvers, such as the quantum Monte-Carlo (QMC) method, provide controlled ap- 
proximations of the single-particle Green's function for lattice models of strongly correlated electrons. 
To access the thermodynamics, however, a thermodynamical potential is needed. We present an effi- 
cient numerical algorithm to compute the grand potential within cluster-embedding approaches that 
are based on novel continuous-time QMC schemes: It is shown that the numerically exact cluster 
grand potential can be obtained from a quantum Wang-Landau technique to reweight the coefficients 
in the expansion of the partition function. The lattice contributions to the grand potential are com- 
puted by a proper infinite summation over Matsubara frequencies. A proof of principle is given by 
applying the VCA to antiferromagnetic (short-range) order in the two-dimensional Hubbard model 
at finite temperatures. 

PACS numbers: 71.10.Fd, 71.27. +a, 71.30.+h 



I. INTRODUCTION 

A powerful strategy to tackle strongly correlated 
fermion systems is to start from a local perspective, i.e. 
to treat the strong local Coulomb interaction exactly for 
an isolated cluster (or atom) first, and to include the cou- 
pling between the clusters in a subsequent step. Dynam- 
ical cluster-embedding approaches^^^ (for reviews 
see Refs. 0Hii3) provide this in a self-consistent way. 
The central object is the cluster single-particle Green's 
function G'(oj) or the cluster self-energy £(w). These are 
computed exactly for the isolated cluster in an external 
dynamical (i.e. frequency dependent) Weiss field. The 
Weiss field mimics the effect of the cluster surrounding 
and is calculated self-consistently using the exact cluster 
quantities and Dyson's equation for the lattice problem 
in the thermodynamic limit. Thereby, one has access 
to the single-particle excitation spectrum, i.e. to the (in- 
verse) photoemission cross section, as well as to different 
physical quantities that can be derived from this. 

The thermodynamics of the system, however, is gov- 
erned by a thermodynamical potential, such as the grand 
potential fl, which is related to the cluster Green's func- 
tion and self-energy in a dynamical cluster embedding 
approach via: 



= fy + TrlnG-TrlnG' 



(1) 



Here fi' is the grand potential of the cluster reference 
system, and G the approximate Green's function of the 
original lattice model in the thermodynamic limit which 
is obtained from the lattice Dyson equation using the 
cluster self-energy. The trace refers to both, spatial and 



temporal lattice degrees of freedom, i.e. involves sums 
over lattice sites and Matsubara frequencies. 

Since dynamical cluster-embedding is a concept that 
directly works in the thermodynamic limit, one of the 
main intentions is to construct phase diagrams. There 
are several well-known situations where the sole knowl- 
edge of the Green's function is insufficient, and a ther- 
modynamcial potential is required. 

A prime example is the Mott transition in the param- 
agnetic phase of the half-filled Hubbard models Using a 
plaquettc of four sites, cellular dynamical mean-field the- 
ory (C-DMFT) predicts a finite £/-range where the metal- 
lic and the Mott insulating solutions are coexisting^ 
While the boundaries of the coexistence region could be 
mapped out precisely by using continuous-time quantum 
Monte-Carlo (CT-QMC)^^ the actual trend of the line 
of first-order transitions has not yet been determined. 
Recent calculations within the variational cluster approx- 
imation (VCA) by using the Lanczos method indicate 
that the first-order line does not end in a second-order 
critical point at zero temperature^ 

Coexistence and competition of phases with different 
long-range or short-range order is also characteristic for 
transition-metal oxides and cuprate materials in partic- 
ular. Salient features of this physics are captured by C- 
DMFT and VCA calculations for the doped single-band 
Hubbard modele d 6 ' 17 Phase separation in the doped 
Mott insulator at T = 0, for example, is another prime 
example where the knowledge of a thermodynamical po- 
tential is necessarily required] 18 i 19 

Note that for cluster-embedding approaches based on 
the self-encrgy-functional theor y 20 : 21 an efficient cvalua- 
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tion of the grand potential via Eq. {1} is decisive not only 
in the final step but during the actual calculation for ex- 
ploiting a variational principle of the form <5f2[E] = to 
find the physical self-energy. Again, the computation of 
f2[S] essentially follows Eq. (fTJ). Furthermore, away from 
the stationary point, metastable phases as well as pre- 
cursors of stable phases can be made visible by looking 
at the functional f2[S] (see Refs. H9II21II22I for examples). 

The cluster size that is accessible using full diagonal- 
ization (ED) or the Lanczos method^ as "cluster solvers" 
is strongly limited. Hubbard clusters consisting of more 
than, say, 12 sites cannot be treated conveniently in this 
way for zero and for finite Therefore, if moderately 
large clusters are needed, finite-temperature quantum 
Monte-Carlo approach (QMC) represents the method of 
choice . 13 i 14 i 25 An important advantage of the stochastic 
method is that uncorrelated ("bath") sites, which make 
up the dynamic Weiss field in the cluster-embedding con- 
text, can be attached to the cluster of correlated sites es- 
sentially without any extra numerical cost. At the same 
time the bath helps to attenuate the QMC sign problem. 

As concerns the evaluation of Eq. ([1]) and thus the 
accessibility of the system's thermodynamics, however, 
ED turns out to be much more convenient: Q' can be 
obtained from the cluster many-body eigenenergies. The 
trace in the third term on the r.h.s. can easily be put in a 
form which involves the poles and the zeros of the cluster 
Green's function G' only, as has been shown in Ref. l2ll . 
The lattice contribution via the trace in the second term 
on the r.h.s. can also be evaluated with this information 
at hand by means of the so-called Q-nmtrix technique.— 

On the other hand, using QMC, there are two main 
obstacles that prevent a straightforward evaluation of 
Eq. {l}: (i) As the Monte-Carlo technique is designed 
to provide expectation values, the cluster grand poten- 
tial fi' cannot be computed directly at finite T because 
of the entropy term. An alternative is to find O' from 
dCl' = -S'dT - N'dfi + D'dU by integrating the cluster 
double occupancies D' over U for fixed temperature T 
and chemical potential fi. This, however, requires sim- 
ulations for a finite U range, (ii) The evaluation of the 
traces must be performed differently since QMC does not 
provide the poles and weights of the one-particle excita- 
tions. Furthermore, the Green's functions G' and G are 
available on the imaginary Matsubara frequencies only. 
Integration along the real axis (as used e.g. in Ref. [22| for 
VCA studies) is thus not possible. Instead, a direct nu- 
mercal summation over the Matsubara frequencies must 
be employed, similar to an integration along the imagi- 
nary frequency axis as described in Ref. [26| for the T = 
limit. 

The purpose of the present paper is to demonstrate 
that these difficulties can be overcome. We suggest to em- 
ploy the continuous-time quantum Monte-Carlo method 
(CT-QMC) 13 i 14 combined with a quantum version of the 
Wang-Landau algorithm i 27 i 28 This allows for a direct nu- 
merical estimate of the cluster partition function for finite 
temperatures. The sign problem remains unaffected. We 



also discuss the evaluation of the traces in Eq. (JTJ by 
summing over Matsubara frequencies. The paper is or- 
ganized as follows: In Sec.|TT]we show how to combine the 
quantum Wang-Landau algorithm with CT-QMC to de- 
termine the grand potential of an isolated cluster. For ac- 
curacy checks, the results are compared with those from 
full diagonalizations of small clusters. In Sec. IIIII the 
Matsubara-frequency summation is discussed emphasiz- 
ing the analytical treatment of the high-frequency limit. 
The application of the technique to the 2D square Hub- 
bard system by embedding a 4 x 4 cluster in an infinite 
square lattice is demonstrated in Sec. IIV1 where various 
thermo dynamical properties are discussed. 

II. QUANTUM WANG-LANDAU APPROACH 

Usual Monte-Carlo algorithms sample configurations 
in the configuration space by an ergodic random walk. 
In practice, however, the Monte-Carlo walk could be 
trapped in a certain part of the configuration space, espe- 
cially when the system is close to a discontinuous phase 
transition. The Wang-Landau algorith m 27 ' 28 has been 
introduced to overcome such problems in the classical 
systems. For quantum systems, it has been proposed 
for an efficient sampling in the context of a stochastic 
series expansion in the inverse temperature^ Basically 
the same idea can be applied to algorithms with different 
expansion parameters. Here we will discuss the appli- 
cation of Wang-Landau algorithm in the context of the 
continuous-time quantum Monte-Carlo metho d 13 i 14 ' 30 i 31 
where the partition function of a fcrmionic quantum sys- 
tem is expanded in powers of the interaction (or hy- 
bridization) strength. 

To discuss this quantum Wang-Landau approach, we 
refer the so-called "weak- coupling" CT-QMC^ 3 . as an ex- 
ample. We also consider the single-band Hubbard model 
which reads 

H = -t ( c L c i<y + c ]o c i°) ~ MX) Hi + U n ^ n H ' 

(2) 

in the usual notation with c irJ being the creation operator 
of an electron at site i and spin projection a =|, J,, with 
ni a = c\ a Ci a7 the nearest-neighbor hopping parameter 
t, the chemical potential fi and the interaction strength 
U. As a non-perturbative and numerically exact method, 
the CT-QMC starts from the infinite sum over diagram 
orders k in the expansion of the partition function: 

~ 00 00 
— =£tfV*) = l + 5>Vfc), (3) 

fc=0 fc=l 

where Zq is the partition function of the non-interacting 
system. The coefficient of the fc-th order 

w(k) = ^w(k,C k ) (4) 

c k 
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is given as a sum over vertex configuration C\ (k vertices 
at order k) which specifies the positions of the vertices 
in space and imaginary time. It includes a fc-dimcnsional 
r-integration over [0, /?]. The weights 

w(k, C k ) = (-Ar/2) fc dct Nl\ k) {C k ) det M (k) {C k ) (5) 

are composed of determinants of matrices Ma k \C k ) 
which are constructed from non-interacting Green's func- 
tions that link the k vertices of a vertex configuration 
Ck- A finite At is used for a discretization of the r- 
integrations but can be made arbitrarily small. This 
sum over vertex configurations is high-dimensional for all 
practical purpo ses . For details of the CT-QMC method 
see Refs. 30, for example. 

A fermion sign problem can be avoided for the model 
at half-filling using a parameter in the interaction 
term, 

H v = — ^[(n iT -a)(n i j.-l+Q!) + (ni T -H-a)(n i |-a)] , 

i 

(6) 

and in the one-particle part 



where q is the momentum vector. It can be further writ- 
ten as 

Z = e W«-«V + e -/9(*,-M+I72)]2 _ (g) 

3 

The high-dimensional sum over orders k and vertex 
configurations C k is sampled by a Monte-Carlo tech- 
nique. Thereby, one can deduce the weight function w(k) 
up to a constant factor. Note that w(k) is not normal- 
ized to unity (Eq. d5J) yields J2kLo w (^) = Z(J=i/Zo)- A 
corresponding histogram, 

p(k) = U k w(k) , (10) 

generated by CT-QMC for a small (2 x 2) Hubbard clus- 
ter, is shown in the left part of Fig. [1] Since w(0) = 
p(0) = 1, as is obvious from Eq. Q, one can in prin- 
ciple determine the unknown factor from the k = 
term in the histogram and find the partition function 
as Z = Zq[1 + YskLiP(k)/p(Q)]' However, this is by no 
means practicable, since usually p(0) is negligibly small 
compared to p(k) at the average order, for example (see 
Fig. [TJ left). Therefore a histogram re- weighting tech- 
nique is necessary. 

As discussed in Ref. the basic idea of Wang-Landau 
sampling is to create a histogram p(k) which is flat for 



all orders k up to a certain cutoff order k c . Thereby, the 
algorithm generates approximately the same number of 
configurations C k both, at low orders and at higher orders 
(up to k c ). To achieve this, a Wang-Landau factor g{k) 
is introduced to re-define the weights w(k, Ck): 

w(k, C k ) -> w(k, C k ) = w(k, C k )/g(k) . (11) 

This also implies the replacement 

w(k) -> w(k) = w(k)/g(k) . (12) 

The (e.g. Metropolis) random walk is performed in the 
usual way, but with respect to the new weights w{k,Ck) 
and thereby with new transition probabilities. The 
Wang-Laudau factor g(k) is chosen to make the new his- 
togram flat, i.e. p(k) = U k w(k) = const for k < k c . 
Hence, the partition function is then given by 

„ oo oo 

— = ^Mfc) = E U k w(k)g(k) . (13) 

fc=0 fc=0 

If the histogram was completely flat to all orders, 

(14) 



(15) 

Note, however, that in practice the Wang-Landau re- 
weighting is performed up to certain order k c only. Fur- 
thermore, in any practical simulation p(k) = const is 
approximate (but can be ensured to arbitrary precision 
in principle). Therefore, Eq. (|13[) with the actual new 
weights and the actual Wang-Landau factor has to be 
used for concrete calculations. The cutoff order k c is basi- 
cally a preselected number representing up to which order 
the histogram will be re-weighted. Typically, k c > (fc), 
where (k) is the average order, has turned out to be a 
good choice. 

As discussed in Refs. [27ll29l . the Wang-Landau factor 
g(k) is continuously modified during a simulation by mul- 
tiplying g(k) with a constant / if a configuration at order 
k is visited, 

g(k) -> g(k) ■ f , (16) 

until the measured histogram p(k) is flat. In practice, 
we set g{k) = e G ^ k \ and increase G(k) at each visit by 
F, i.e. G(k) — > G(k) + F. In our implementation of the 
algorithm, the flatness of p(k) is defined by requiring the 
difference of the smallest and largest p(k) to be within a 
given small tolerance ry > 0, i.e. 

Prain 

> (1 - n)p 

max ■ 

(17) 



— =p(0)X>(*) 



fe=0 



H a = -t E (cLc Jff + C ] ff c, CT )-(M-C//2)E^-^(«-« 2 )^- 

(7) 

Here N is the total number of sites in the system. The and with p ^ = x wc would get 
non-interacting partition function can be determined eas- 
ily by diagonalizing Hq. This yields: 



Z = e /5JV(«-c 2 )c/ Xr £ -(3 -M+i'/a)™,, j rg\ 



Z/Z = J2g(k)/g(0) 



k=0 
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In the actual simulation, multiple steps of Wang-Landau 
re- weighting is recommended. Each time when the above 
condition is fulfilled, a further decrease of the multiplica- 
tive constant / improves the flatness of the histogram. 

The Wang-Landau factor g(k) is positive by construc- 
tion. This implies that the new weights C&) are 
positive if and only if the original weights w(k,C k ) are 
positive. Therefore, the re-weighting technique does not 
introduce a new source for a sign problem. 

Fig. [T] (right) shows the rc-weightcd histogram. We 
can see that after the Wang-Landau re-weighting up to 
k c = 40, the histogram p(k) is sufficiently "flat". Each 
time when the criterion Eq. (fT7|) is fulfilled, F has been 
reduced by a factor two, and the tolerance decreased by 
0.05 to improve the flatness of the histogram. For the 
example shown and also for all calculations below F — 
0.01 and a tolerance r] — 0.2 as final values have turned 
out to be sufficient. Note that if Z/Zq is computed from 
Eq. p3p it is unnecessary to have a strictly flat histogram. 

In Fig. [21 the grand potential £1 = — ThxZ is shown 
for three different cluster sizes. f2 is calculated as a func- 
tion of the chemical potential at (3t = 10 and U/t = 4 
and compared to numerically exact results from full di- 
agonalization of the problem. The agreement with the 
exact grand potential is excellent. The standard devia- 
tion of the cluster grand potential is reasonably small in 
our calculation, which is crucial since the lattice grand 
potential consists of the cluster grand potential f2' and 
of Tr In G — Tr In G' . Both contributions are of the same 
order of magnitude. Hence an accurate determination 
the cluster grand potential is important. Within the 
Wang-Landau CT-QMC implementation, we find that 
the statistical error of Q, 1 is comparable to that of the 
Tr In G Tr In G', and typically three orders of magnitude 
lower than its mean value for the calculations presented 
here. 
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FIG. 1: Left: Histogram p(k) of the order k in the weak- 
coupling expansion. Calculation for a 2 x 2 Hubbard cluster 
at half-filling, j3t = 5 and U/t = 4. Right: Histogram with 
Wang-Landau re- weighting, p(k), choosing k c = 40. 



Besides the possibility to compute the grand poten- 
tial, there is another important advantage of the quan- 
tum Wang-Landau method: Suppose that we have the 
Wang-Landau factor g(k) from a calculation at interac- 
tion strength Uq. Then the full U dependence of the 
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FIG. 2: The grand potential per site as a function of the 
chemical potential /j, for small Hubbard clusters at fit = 10 
and U/t = 4. The CT-QMC results are shown as dots. The 
statistical error is smaller than the symbol size. Full ED re- 
sults are shown as lines for comparison. 



average of an observable can be determined at once for 
< U < Uo since the U dependence of the weights is 
trivially U k w(k,C k ). Since w(k,Ck) is a functional of 
the free Green's function and (3 only, it is U indepen- 
dent. 

The average of an observable O is 

_ E k E Ck u k Hk,c k )Q(k,c k ) 
{0) - E k Ec k u^(k,c k ) ' (18) 

Using Eq. ([Tr]) , we have 



(O) 



J2 k u k g(k)J2 Ck w(k,c k )0(k,c k ) 



Y Jk U k g{k)Y jCk w{k 1 C k ) 
Define an average at a given order k by 

j: Ck w(k,c k )o(k,c k ) 



0(k) 



Ec k ™(k,C k ) 



(19) 



(20) 



This average is carried out by importance sampling of 
configurations C k in the Monte-Carlo technique. The 
subsequent average over different perturbation orders k 
then reads as 



(O) 



E fc U k g{k)Q{k) 
E k U k 9(k) 



(21) 



The U dependence of (O) is now explicit. 

The [/-dependent average perturbation order (fc), the 
average particle numbers, the interacting Green's func- 
tion, etc. are easily obtained from this equation. Note 
that since U/Uq < 1, only smaller number of diagrams is 
required in the construction of the full partition function 
at U compared to Uq, while g{k) is determined up to the 
maximum contributing order at Uq- 
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Fig. [3] shows the Wang-Landau factor g(k) at fit = 5 
and U/t = 4 on a 2x2 cluster as well as the on-site 
Green's functions for U — 0.5 to U = 4 with steps 
AU = 0.5 as functions of the Matsubara frequency. The 
Green's functions have been calculated using Eq. (|21[) . 
We have checked that these results agree with those from 
conventional CT-QMC simulations at the respective fixed 
U. The low-frequency behavior is reminiscent of the 
metal-insulator transition in the infinite square lattice. 




FIG. 3: The on-site Green's function for a 2 x 2 cluster at fit = 
5. From bottom to top, different lines corresponds to different 
interaction strengths varying from U/t = 0.5 to U/t = 4.0 in 
steps of 0.5. The Wang-Landau factor g(k) was determined 
at U/t — 4, and the Green's function for other interaction 
strengths are calculated from Eq. (12111 , 



III. 



LATTICE CONTRIBUTION TO THE 
GRAND POTENTIAL 



Within a cluster-embedding approach, the approxi- 
mate grand potential of the lattice fermion model in the 
thermodynamical limit with parameters t and U is given 
via Eq. {T]) by the grand potential f2' of a small cluster 
with parameters t' and U and a lattice contribution. The 
latter is obtained from the approximate single-particle 
lattice Green's function, G^\j = iu> n + /i — t — 5V,u", and 
the exact cluster Green's function, G^,\j = iuj n +/j, — t' — 
E t >u, as 



Tr]n(Gt, u )-T[]n(G t ,,u) 



2 r y e -,o+^y trln ^ 



G t >,u(k,iwT, 



N 



G' tlu (iu) n ) 



N ^ 



e l "" 0+ lndet[l - V~ k G v ,u(^n)\ (22) 



Here Tr = 2T^g e IW " 0+ tr combines the usual trace 
referring to the cluster sites, the sum over Matsubara fre- 
quencies, and the sum over wavevectors in the reduced 



Brillouin zone of the superlattice of disconnected clus- 
ters. The factor 2 accounts for the two spin projections. 



Furthermore, Vt 



t! is the hybridization matrix 



stressing the difference between the lattice and the clus- 
ter system with respect to the one-particle parameters. 
N c and N are the size of the cluster and of the lattice, 
respectively. 

There are in principle two different techniques that can 
be employed to evaluate the frequency sum in Eq. (j22j) : 
(i) An analytical evaluation is possible, once the poles 
and the weights of the cluster Green's function are avail- 
able. This is the idea of the Q-matrix metho d 21 ' 32 which 
works well at zero temperature where only a few num- 
ber of poles contribute, and for small clusters where the 
full eigensystem is available from exact diagonalization. 
It is obvious, however, that the CT-QMC technique is 
not capable to locate the poles, (ii) Alternatively, one 
can evaluate the frequency summation numerically us- 
ing the CT-QMC result for Gt> ui^n) on the Matsubara 
frequencies. This summation can be performed in close 
analogy to the integration along the imaginary frequency 
axis described by Scncchalj^ 

Below we briefly describe the method that has been 
used here: Note that a direct numerical frequency sum- 
mation, in Eq. (|22p . is not possible since the factor e luJn0+ 
is crucial for convergence. One therefore has to treat 
the high-frequency part separately and analytically. In- 
troducing a (sufficiently large) cutoff frequency lua, the 
infinite Matsubara sum can be split up in the following 
way: 

H\a{G, u )-H\a{G,, u ) 



N 

k n=-A 



3 *c„o+ ln det [l-V k G t >,u(iu n )] 



det [l - Vk/iu>n\ 



- 2T £^ E 



N 

k n=-oo 



In det 



(23) 



The first term involves a finite Matsubara sum only and 
can thus be computed by direct numerical summation 



with 



set to unity. For the second term, the fre- 



quency summation can be done analytically. We find: 



-2TY^ Y e -" 0+ lndct 

k n— — oG 



Vi 



= 2T 



N 



(24) 



where F- aa are the diagonal elements of Vjr and a refers 
to the sites in the cluster. This yields 

^jV c A detTl- VtG t >u(iij n )] 
fi = fi'-2TV— V In- L 

N 

k n=-A 

+2TN c ln2-2Tj2 — ln(l + e- 9 i"/ T ) (25) 



det [l — Vj,/icu n ] 



a . k 
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for the lattice grand potential of a cluster-embedding 
approach at finite temperatures. Within CT-QMC, the 
cluster contribution 1Y is determined as shown in Scc.HH 
Note that, for different cluster-embedding methods, dif- 
ferent implementations of the non-interacting partition 
function Z must be considered (see Ref. HH) . 

The one-particle Green's function Gt',u(iw n ) is mea- 
sured up to the cutoff frequency u>\. Clearly, the accu- 
racy of the frequency summation in the lattice contri- 
bution to the grand potential is sensitive to this cutoff. 
To check the accuracy, we consider as a simple test case 
the half-filled Hubbard model with a semi-elliptical free 
density of states po(z) of bandwidth W = 4i, and the 
Hubbard atom (N c = 1) without bath sites as a reference 
"cluster". The cluster Green's function at half-filling is 
readily obtained, the corresponding local self-energy is 
5V t U (iui n ) — U 2 1 A%Lo n . Within the cluster-embedding 
method, this approximates the self-energy of the lattice 
model and yields the lattice Green's function, the local 
elements of which are thus given in terms of the free den- 
sity of states po(z) as 



W/2 

-W/2 lu} n 



po(z)dz 



-'t'U 



(26) 



The poles of lattice Green's function are determined by 
solving the equation ioj n — z — Tf ,u{i^n) = 0. From 
these and from the poles of G t > u(iu> n ) the lattice grand 
potential can be determined analytically. We find: 



il = SI' + T In 



1 + cosh( 



2T' 



T / p (z)dz\n 



. z , Vz 2 + U 2 
cosh — + cosh — 



(27) 



This is easily evaluated numerically to any desired accu- 
racy by means of the Simpson method. Table U compares 
the exact result for the lattice grand potential (per site) 
with the numerical result obtained by evaluation of Eq. 
(f25|) for different cutoff frequencies. This demonstrates 
that a moderate cutoff frequency is sufficient to obtain a 
reliable result with a relative error, as compared to the 
exact result from Eq. (f2"7]) . of the order of 10 -3 . This 
improves with increasing cutoff frequency. With decreas- 
ing temperature, the required frequency cutoff becomes 
larger. 





Exact 


Numerical: Eq. ([251 


1.0 


-0.03423 


A = 10 


A = 30 


-0.03425 (5.8E-4) 


-0.03424 (3.0E-4) 


10.0 


-0.03099 


A = 50 


A = 100 


-0.03106 (2.3E-3) 


-0.03102 (9.7E-4) 



TABLE I: Grand potential per site for the Hubbard model 
at half-filling with semi-elliptical free density of states (band- 
width W = 4i) as obtained for a reference "cluster" with 
N c = 1. Numerical results obtained from Eq. (|25[1 for differ- 
ent temperatures and cutoff frequencies are compared to the 
exact result (Eq. (|27[lh The relative differences are given in 
brackets. 



cluster-embedding approach and thus sufficient to pro- 
vide a proof of principle. Clearly, our method to deter- 
mine a thermodynamical potential can easily be extended 
to study larger clusters and lower T with more computa- 
tional effort. 

Generally, the additional computational effort neces- 
sary for the rewcighting technique is reasonable: First, 
we note that the time needed to initially determine the 
Wang-Landau factor is usually less than the time spent 
for the measurement of obscrvablcs in our calculations. It 
is even negligible in those cases where there is a good ini- 
tial estimate for the Wang-Landau factor, as e.g. from 
a preceeding calculation with slightly different model 
parameters. Second, since the reweighting technique 
amounts to importance sampling for any given order k 
but full sampling over all k up to the cutoff order, the 
configuration space is enlarged considerably. However, 
the orders to be sampled additionally arc computation- 
ally inexpensive because they are significantly lower than 
the average order. The numerical effort is still dictated 
by the average order as in conventional CT-QMC. 

To keep things simple we furthermore use clusters 
without bath sites, i.e. we perform a VCA calculation 
where a single variational parameter is considered only. 
Note that bath sites could be added without any ad- 
ditional numerical cost within the weak-coupling CT- 
QMC. For the particle-hole symmetric model at half- 
filling, there is no sign problem. 

The strength h of a staggered magnetic Weiss field sug- 
gests itself as the most relevant variational parameter. 
Hence, we add the term 



Weiss 



(28) 



IV. SHORT-RANGE ANTIFERROMAGNETIC 
ORDER AT FINITE TEMPERATURES 

For an application of CT-QMC combined with the 
quantum Wang-Landau method we consider the Hub- 
bard model on the square lattice at half-filling and j3t < 5 
and use a 4 x 4 Hubbard cluster for the embedding. This 
is well beyond the cluster size that can be accessed (con- 
veniently) by exact diagonalization or Lanczos within a 



to the Hamiltonian of the reference cluster. Here, Q = 
(■7T, 7r) is the antifcrromagnetic wave vector. To get the 
optimal value of the Weiss field h within the framework 
of the self-energy-functional approach, we compute the 
lattice grand potential f2 for different h according to the 
method outlined in the previous sections and search for 
stationary points: 

dCl 
~dh 



. 



(29) 
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FIG. 4: The internal energy per site determined from a 4 x 4 
embedded cluster at different temperatures. From bottom 
to top, the interaction strength increases from U/t = 1.0 to 
4.0. Statistical error bars are smaller than the symbol size 
and hardly visible. Our results are in good agreement with 
the QMC solution on an 8 x 8 cluster obtained in Ref. [34] at 
U/t = 4.0. 



This optimization of the thermodynamical potential can 
also be performed in the presence of a finite physical ex- 
ternal staggered field of strength /iaf which must not be 
confused with the variational paramerter h. The AF or- 
der parameter, i.e. the staggered magnetization, is then 
obtained as 



ssi 

llm "HI — 

h A F^0 0/1 AF 



(30) 



Spontaneous symmetry breaking is indicated by a finite 
to, or, looking at Sl(/i) for hxv = 0, by a finite optimal 
value for the variational parameter h. 




FIG. 5: The grand potential Q, per site from Eq. (|25[) and the 
cluster grand potential Q' per site as functions of temperature 
for two different interaction strengths. 

For our calculations we usually re-weighted the his- 



tograms up to the order k c — N c [3U/2. For higher tem- 
peratures, it turned out to be useful to extend k c by 20 
more orders for the re-weighting. The Wang-Landau fac- 
tor g(k) was kept fixed when the histogram for k < k c 
became sufficiently flat, as controlled by Eq. (fl"7[) . 

Each observable was measured by totally 2.56 x 10 8 
Monte Carlo steps after a flat histogram has been ob- 
tained, see Eq. (jTTJ) . This turned out to be sufficient 
for controlling the statistical error. In most of the re- 
sults shown below, error bars are smaller than the sym- 
bol size. The frequency summations discussed in Sec. IIIII 
have been carried out with a frequency cutoff of A = 120 
for all temperatures considered here. 



1.4 
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0.8 
0.6 
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FIG. 6: Entropy (per site) as calculated from Eq. (|32[) for the 
same parameters as for Fig. [4] and [5] 



To start with, we discuss results for vanishing Weiss 
field, i.e. h = 0. Fig. Q] shows the internal energy 
as a function of temperature for different interaction 
strengths. The internal energy is computed as E = 
Ekin + U(D) from the double occupancy (D) and the 
kinetic energy Ekin which is obtained from the lattice 
Green's function (i.e. with self-energy replaced by the 
cluster self-energy) <i The agreement of our results at 
U/t = 4.0 with those of a QMC calculation^! for an iso- 
lated Hubbard cluster of 8 x 8 sites is excellent in the 
entire temperature range. 

While the internal energy (as an expectation value of 
an observable) is directly available within usual QMC 
methods, a thermodynamical potential is not. The ad- 
vantage of the Wang-Landau approach lies in the direct 
accessibility of the partition function and thus of the 
grand potential for any U < Uo — if implemented within 
the (weak-coupling) CT-QMC framework. 

Fig. shows the temperature dependence of the grand 
potential S7 as obtained from Eq. ([23)) for two differ- 
ent interaction strengths. The cluster grand potential 
SI' = —T\nZ' with Z' determined by the Wang-Landau 
approach is shown in addition. Using the 4x4 cluster, 
SI' already represents the main contribution to the to- 
tal grand potential. As a function of temperature, the 
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lattice contribution Tr In G — Tr In G' resulting from the 
cluster embedding is basically constant at high tempera- 
tures but is of crucial importance for low T. This has to 
be expected as spatial correlations grow with decreasing 
T. 




from Eq. (32) 
differentiation of £1 



2.5 
T/t 



FIG. 7: The comparison of the entropy from Eq. (|32|l and the 
numerical differentiation of the the lattice grand potential, 
S — —dQ/dT, at U/t = 1 for different temperatures. 

The general trend of f2(T) is ruled by general thermo- 
dynamical relations. Its negative slope is given by the 
entropy, the second derivative corresponds to the specific 
heat divided by T . The linear high-temperature trend 
reflects the Hilbert-space dimension A N of the Hubbard 
model with N sites since 



lim S(T) = ln4 N = 27V In 2 

T — *oo 



(31) 



for the entropy. At low temperatures, the results might 
be indicative of a quadratic behavior around T = 0, cor- 
responding to a linear entropy and a linear specific heat, 
as one would expect for a Fermi liquid. Note that there 
is hardly a change of the temperature trend when vary- 
ing U. More definite statements, however, would require 
calculations at considerably lower T. 

Fig. [S] shows the entropy calculated from 



E-n- fiNn 
T 



(32) 



as a function of temperature for two different interaction 
strengths. Here n is the particle number per site. It 
is obvious that at low T and for U/t — 4 the results 
for S are less accurate when compared with those for the 
grand potential or internal energy although the statistical 
error is reasonably small. Via Eq. I|32p. a small T in the 
denominator enhances the systematic error of E and f2 
resulting from the cutoff k c . To get a more satisfactory 
result, k c would have be to increased to include more 
diagrams in the CT-QMC. 

The entropy can be also obtained as S = —dVl/dT 
from the grand potential. There is no reason why the two 



ways for computing S should give the same result within 
an approximate theory. Opposed to quantities like the 
particle number or the staggered magnetization (see Ref. 
|35|. for example), there is no internal consistency of a 
cluster-embedding approach with respect to the entropy. 
We have therefore explicitly checked for differences be- 
tween both. Fig. [7] shows the temperature dependence of 
the entropy at U/t = 1 as obtained from Eq. (f3"2")) and 
from Fig. [5] by differentiation, respectively. Finite differ- 
ences show up at intermediate temperatures. To a major 
extent, these are due to the error in calculating the tem- 
perature derivative from the discrete and small number 
of O values. This can be seen e.g. at temperature T/t = 2 
in the figure. There is, however, a small but significant 
remaining difference between both ways for computing 
S. This represents an intrinsic problem of the theory, 
and actually of any cluster-embedding approach. The 
observed inconsistency may serve as a measure for the 
error due to the finite cluster size since thermodynami- 
cally consistent results for S can be expected strictly in 
the infinite-cluster limit only. 

Finally, Fig. [H] demonstrates that our technique can in 
fact be used for the determination of variational param- 
eters within the framework of the self-energy-functional 
theory. The figure shows the lattice grand potential f2 as 
a function of the strength of the Weiss field h (see Eq. |2"8")) 
for different U and T. Let us concentrate on the results 
obtained for the 4 x 4 cluster first. For U/t = 1, U/t = 2 
and fit = 5.0 (panel a) we find a single stationary point 
at h = 0. This is indicative of the paramagnetic phase 
at high temperatures and weak interaction. For U/t = 4 
(panel a) , the SFT grand potential clearly displays a min- 
imum around h = 0.15. This corresponds to antifcrro- 
magnetic order. Here we also get a non-zero value for the 
order parameter from Eq. (|3"0)) . Note that the variation 
of il with h is small and comparable to the statistical 
error. This shows that (3t = 5 is close to the Neel tem- 
perature for U/t = 4 (and for the given cluster size, see 
below). As the trend of fl(h) is symmetric with respect 
to h, the point h = represents another stationary point 
corresponding to the paramagnetic phase. The latter is 
metastablc as Vt(h = 0) is higher than 0(/i = 0.15). 

Fig. [8l panel (b) displays results for a higher temper- 
ature (f3t = 2). Here we are again left with the param- 
agnetic phase for all U/t only. Obviously, the variation 
of the grand potential with h is most pronounced for 
U/t = 4 while is becomes more and more flat with de- 
creasing interaction strength. This is due to the fact that 
the h dependence enters the self-energy functional via the 
self-energy only, i.e. ft(h) = 0[S(/i)], and fi[X] = in the 
non-interacting limit. For fit = 5 (panel a), this is dif- 
ferent. Here, the trend of fl(h) first becomes stronger 
(comparing U/t = 1 with U/t = 2) as explained above. 
For stronger interactions, however, this mechanism has 
to compete with the tendency to form a minimum at a 
finite h. This explains the non-monotonic trend in U 
visible in panel a. 

Due to the Mermin- Wagner theorem^ there is no long- 
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FIG. 8: Lattice grand potential f2 as a function of the variational parameter h as obtained by embedding a 4 x 4 cluster (a,b) 
and a 2 x 2 cluster (c). Results for different interaction strengths U/t as indicated in the figure. Note that the difference 
tt(h) - fi(0) is plotted, (a) fit = 5.0, N c = 4 x 4, (b) fit = 2.0, iV c = 4 x 4, (c) fit = 5.0, N c = 2 x 2. Symbols with statistical 
errors represent the numerical data. Lines are obtained by a spline interpolation and serve as a guide to the eye only in panel 
(a) and (b). Lines in panel (c) are from full ED calculations for comparison. 



range antiferromagnetic order in the two-dimensional 
Hubbard model at finite temperatures. Therefore, in 
a strict interpretation, the symmetry-broken phase ob- 
tained from the cluster-embedding approach has to be 
seen as a mean- field artifact.. In the context of the 
dynamical cluster approximation, for example, this has 
been studied extensively in Ref. 1371 . For our case, a non- 
zero value of optimal Weiss field actually indicates that 
the antiferromagnetic correlation length £ exceeds the lin- 
ear size of the cluster, i.e. four sites. 

This interpretation is corroborated by the comparison 
of the results from the 4x4 cluster embedding with the 
corresponding ones obtained from a 2 x 2 cluster (see 
panel c). Here we have also added the results using full 
ED as a cluster solver for comparison. Antiferromag- 
netic (short-range) order at finite temperatures is build 
up with increasing interaction. If, within a cluster mean- 
field approach, a symmetry-broken state is obtained once 
£ ~ V^Vc, one should therefore expect a lower critical 
interaction strength when reducing the cluster size iV c . 
This is exactly what is found when comparing the re- 
sults for N c = 4 x 4 (panel a) with those for iV c = 2 x 2 
(panel c). At U/t = 4 the optimal value for h for for 
the 2x2 cluster (h opt = 0.24) is larger than that for the 
4x4 cluster (h opt = 0.15). Furthermore, the difference 
fl(h) — f2(0) , which measures the stability of the antifer- 
romagnetic order, is higher for the 2x2 calculation as 
compared to the 4x4 one at the same U/t. 



V. CONCLUSIONS 

Besides the calculation of one-particle excitations, 
cluster-embedding approaches are able to provide de- 
tained information on the thermodynamical properties of 
correlated lattice-fermion models, such as the Hubbard 
model. When combined with the quantum Monte-Carlo 
technique as a cluster solver, however, cluster DMFT 
schemes or the VCA cannot directly access a thermo- 



dynamical potential because at finite temperatures the 
entropy term in the free energy, for example, cannot be 
obtained by measuring an observable. On the other hand, 
a thermodynamical potential is of crucial importance for 
the construction of phase diagrams as it decides on the 
relative stability in situations where there are different 
competing phases, e.g. in the high-T c problem. 

For the novel continuous-time QMC schemes, the 
present study has demonstrated that there is an elegant 
way to overcome this difficulty, namely by combining the 
CT-QMC approach with a quantum Wang-Landau tech- 
nique. Employing a proper reweighting of the transition 
probabilities, it is easily possible to construct a flat dis- 
tribution of the perturbation order k up to an in principle 
arbitrarily high cutoff order. This allows for a direct cal- 
culation of the finite- T partition function of the cluster. 
This has been demonstrated here for the weak-coupling 
variant of CT-QMC but can be generalized to other 
schemes. Note that modest changes of the QMC code are 
requricd only. Another advantage of the Wang-Landau 
technique which is worth mentioning consists in the fact 
that the U dependency of obscrvables is obtained in a 
single simulation: Since the weight factors are available 
for any perturbation order separately, the U -dependence 
of obscrvables is trivially given within the weak-coupling 
CT-QMC. Once the perturbation-order cutoff k c is suit- 
ably fixed for a calculation at interaction strength £/o, 
the same weight factors determine the entire U depen- 
dence of the observables for U < IV Note, however, 
that one cannot make use of this advantage within (clus- 
ter) DMFT, since the effective cluster problem is itself U 
dependent. 

The lattice contribution to the total grand potential is 
obtained from the cluster Green's function and the lat- 
tice Dyson equation by a proper summation over Mat- 
subara frequencies in the TrlnG-likc terms. To obtain 
convergent results, the high-frequency asymptotics must 
be controlled carefully. 

Compared to full diagonalization of the cluster prob- 
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lem or to Lanczos-type approaches, the CT-QMC- Wang- 
Landau method can used for larger clusters with mod- 
erate computational effort. For a proof of principle, we 
have applied the VCA to the Hubbard model on a square 
lattice at half-filling and considered a single variational 
parameter only. A cluster with 4x4 sites has been used 
for the embedding. This is clearly beyond the range that 
has been accessible in former VCA studies (at T = 0) 
which were based on the Lanczos method and thereby 
restricted to clusters with about 12 sites at most. 

It is important to note, however, that one of the ad- 
vantages of the QMC technique is that uncorrelated bath 
sites can be integrated out and can thus be included in 
the cluster reference system without any extra numeri- 
cal cost. The calculations presented here have been done 
without bath sites. While this is sufficient to discuss the 
methodical aspects of the CT-QMC- Wang-Landau ap- 
proach, future studies should be carried out by including 
a continuous bath, i.e. the VCA should be replaced by 



cellular DMFT (or a different cluster-DMFT variant) for 
finite-T calculations based on QMC. This ensures an op- 
timal description of the local quantum fluctuations on the 
cluster without extra cost. The implementation of the 
quantum Wang-Landau technique as well as the compu- 
tation of the lattice contribution is identical to the VCA- 
based study shown here. 
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